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Preface 
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aspects  of  the  project,  my  wife  Marisa  was  a  true  inspira¬ 
tion.  Her  thoughtful  support  was  priceless  to  me.  And  to 
Ms.  Phyllis  Reynolds  my  very  capable  typist,  I  owe  my  most 
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I  Abstract 

V 

Four  numerical  methods  are  used  to  solve  a  spe¬ 
cific  set  of  problems  and  then  the  methods  are  compared  for 
accuracy  and  efficiency. 

The  methods  are:  standard  finite  differences, 
collocation,  Galerkin,  and  least  squares.  The  latter  three 
methods  are  finite  element  methods  which  use  either  Lagrange 
linear,  Hermite  cubic,  or  Hermite  septic  piecewise  poly¬ 
nomials  as  interpolation  functions. 

The  problem  set  consists  of  second-  and  fourth- 

v.  ^ — ''' 

order,  linear  and  nonlinear,  differential  equations  with 
constant  and  variable  coefficients.  The  linear  equations 
govern  elementary  structural  members  and  the  nonlinear  equa¬ 
tion  is  a  one-dimensional  analog  for  transonic  flow  past  an 
airfoil. 

The  three  major  conclusions  are:  (1)  the  least 
squares  method  with  Hermite  cubic  polynomials  was  the  method 
of  choice  for  the  second-order  linear  equations,  (2)  the 
collocation  method  was  chosen  over  the  Galerkin  and  the 
finite  difference  methods  for  the  fourth-order  equations, 
and  (3)  Galerkin  method  was  chosen  over  the  collocation 
method  for  the  nonlinear  problem. 

A 
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COMPARISON  OF  NUMERICAL  ANALYSIS  METHODS  FOR 


SOLVING  ONE-DIMENSIONAL,  ELLIPTIC 
DIFFERENTIAL  EQUATIONS 


I ,  Introduction 


Background 

Comparisons  of  numerical  analysis  methods  for 
solving  differential  equations  have  been  of  interest  for 
many  years.  One  of  the  most  recent  studies  was  performed 
by  Houstis  et  al.  (Ref  5) .  In  this  study  four  methods  were 
evaluated  for  solving  second-order,  linear,  elliptic, 
partial  differential  equations.  The  four  methods  were: 
standard  finite  differences;  collocation;  Galerkin;  and 
least  squares  using  Hermite  cubic  piecewise  polynomials. 

They  concluded  that: 

1.  There  is  normally  a  "crossover  point"  at  low  accu¬ 
racy  beyond  which  collocation  is  more  efficient  than 
standard  finite  differences.  Even  when  finite 
differences  is  more  efficient,  it  is  by  a  small 
amount  while  collocation  is  sometimes  dramatically 
more  efficient  than  finite  differences.  Colloca¬ 
tion  is  much  superior  for  problems  whose  boundary 
conditions  involved  derivatives. 

2.  There  is  practically  no  difference  at  all  between 
Galerkin  and  least  squares  performance.  They  tend 
to  be  slightly  more  accurate  than  collocation  but 
are  very  much  less  efficient  because  of  the  increased 
work  to  compute  the  coefficients  in  the  matrix  prob¬ 
lem  to  be  solved  [Ref  5:324-325]. 


In  view  of  this  work  the  question  arises,  will  these  results 
hold  for  higher-order  linear  and  nonlinear  differential 
equations? 

The  three  differential  equations  selected  for  com¬ 
parison  analysis  in  this  thesis  are  equations  which  govern 
elementary  problems  from  the  field  of  aeronautical  engineer¬ 
ing.  The  first  problem,  the  axial  deflection  of  a  rod,  is 
governed  by  a  second-order,  linear,  ordinary  differential 
equation.  The  second  problem,  the  bending  of  a  be.n,  is 
governed  by  a  fourth-order,  linear,  ordinary  differential 
equation.  In  practice  these  two  structural  members  are 
used  to  model  more  complicated  structures.  For  example, 
in  the  computer  program  "ANALYZE,"  used  by  the  Air  Force 
Flight  Dynamics  Laboratory  for  analysis  of  aerospace  struc¬ 
tures,  the  rod  member  is  used  to  model  spar  and  rib  caps 
and  other  line  elements  (Ref  11:18).  The  third  differential 
equation  is  second-order  and  nonlinear;  it  represents  a 
one-dimensional,  steady,  analog  to  the  problem  of  transonic 
flow  over  an  airfoil.  Fung,  et  al.  (Ref  3),  solved  this 
equation  with  unsteady  perturbations  using  a  finite  differ¬ 
ence  method.  The  finite  element  methods  included  in  this 
thesis  have  not  been  used  previously  to  solve  this  non¬ 
linear  equation.  The  purpose  of  this  thesis  is  to  determine 
if  the  results  of  Houstis  (Ref  5)  apply  to  these  three  prob¬ 
lems. 

The  numerical  methods  evaluated  by  Houstis  and  used 
in  this  study  can  be  classified  into  two  categories:  finite 
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difference  methods,  and  finite  element  methods.  Three  formu¬ 
lation  techniques  were  used  to  obtain  the  finite  element 
equations.  They  are:  collocation,  Galerkin,  and  least 
squares  methods.  These  numerical  methods  as  used  in  this 
study  transform  the  original  continuous  differential  equa¬ 
tion  into  a  system  of  algebraic  equations.  These  equations 
can  be  solved  to  yield  an  approximate  solution  to  the 
original  continuous  problem.  Although  the  four  methods  all 
yield  algebraic  equations,  the  approaches  are  very  different 
for  finite  element  methods  as  compared  to  finite  difference 
methods. 

The  finite  difference  method  approximates  the 
derivatives  appearing  in  the  governing  differential  equa¬ 
tion  by  difference  quotients  (Ref  4:222-228).  These  quo¬ 
tients  are  formulated  at  N  (a  finite  number)  points  in  the 
problem  domain.  This  process  generates  N  equations,  non¬ 
linear  in  general,  with  N  unknowns.  After  implementing  the 
boundary  conditions,  the  resulting  equations  can  be  solved 
for  the  remaining  unknown  values.  These  values  of  the 
function  at  N  points  form  a  discrete  approximation  to  the 
solution  of  the  problem. 

Unlike  the  finite  difference  method,  the  finite  ele¬ 
ment  methods  yield  a  piecewise  continuous  approximation  to 
the  solution  (Ref  6:3-9).  The  approximation  is  achieved 
by  discretizing  the  domain  into  E  elements,  with  N  nodes 
and  M  nodal  parameters.  The  approximate  solution  in  a 
"global"  sense  is  given  by: 
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(1) 


4>  = 


M 

E  Vi 

j-=i 


where 

<f>  =  global  approximate  solution, 

Nj  -  global  shape  (interpolation)  functions,  and 
<t>j  =  nodal  parameters  (field  variables) 

When  Eq  (1)  is  substituted  into  the  continuous  differential 
equation  there  will  be  some  error  for  all  cases  except 
where  the  number  of  parameters  M  is  infinite,  or  where  the 
shape  functions  contain  all  terms  of  the  exact  solution. 

The  error  e  is  given  by 

e  =  £(*)  -  f  (2) 


where 

e  =  vector  of  nodal  errors, 

£  =  differential  operator,  and 
f  =  vector  of  equivalent  nodal  forces. 

The  finite  element  equations  are  derived  in  this  thesis  by 
the  method  of  weighted  residuals  (Ref  6) .  This  method 
requires  the  residual  error  e  be  zero  in  some  average  sense. 
The  generalized  orthogonality  condition  is  used  to  achieve 
this  result.  In  general,  the  error  is  forced  to  satisfy  the 
following  equations 

/ei^dfi  =  0  i=l,2, . . .  ,M  (3) 

u 
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where 


ft  =  problem  domain, 
ij>^  =  weighting  functions. 

The  procedure  yields  M  algebraic  equations  with  M  unknown 
nodal  parameters.  When  solved  for,  the  parameters  can  be 
substituted  in  Eq  (1)  to  give  the  approximate  solution  to 
the  problem.  These  equations  can  be  simplified  by  writing 
Eq  (3)  as  the  sum  of  integrals  over  each  element.  By  using 
an  appropriate  coordinate  transformation  each  integral  can 
be  evaluated  on  a  local  basis.  Likewise,  the  approximate 
global  solution  4>  is  assumed  to  be  the  sum  of  the  approxi¬ 
mate  solutions  within  each  element. 

Thus 

♦  -  y,  *<e)  <4) 

e=l 

(e) 

where  <J)  '  is  the  approximate  "local"  solution  in  element 
e.  These  local  solutions  are  related  to  the  elemental 
nodal  parameters  by  the  following  equation: 


m 


(e) 


=  S  N.c 
j=l  3 


(5) 


where 


m  =  number  of  nodal  unknowns  associated  with 
element  e, 

N.  =  local  shape  functions  defined  for  element  e, 
1  and 


=  nodal  parameters  associated  with  element  e. 


Then  if  the  residual  error  in  each  element  is  forced  to 
satisfy  Eq  (3)  over  the  element's  domain,  the  result  is  m 
equations  in  m  unknowns.  These  m  equations  will  have  the 
same  form  for  all  like  elements.  Thus,  the  equations  for 
other  like  elements  are  easily  generated,  and  the  global 
equations  can  be  assembled  from  the  E  sets  of  local  equa¬ 
tions.  This  similarity  between  like  elements  greatly 
reduces  the  amount  of  work  required  to  formulate  the  global 
equations.  More  of  the  details  for  the  finite  element  tech¬ 
niques  are  given  in  Appendices  C  and  D. 

Approach 

The  general  approach  followed  by  Houstis  (Ref  5) 
is  used  in  this  work.  This  approach  consists  of  choosing: 
a  problem  set  from  the  specified  domain,  the  numerical 
methods,  and  the  performance  criteria  for  evaluating  the 
methods.  The  next  step  is  to  select  and  solve  specific 
problems;  and,  lastly,  evaluate  each  method  according  to  the 
criteria  established. 

The  problem  set  has  been  briefly  described  in  the 
previous  section.  It  consists  of  second-  and  fourth-order, 
linear  differential  equations  with  both  constant  and  vari¬ 
able  coefficients.  The  last  problem  is  a  second-order 
nonlinear  differential  equation.  The  two  constant  coef¬ 
ficient,  linear  equations  are  solved  for  two  sets  of 
boundary  conditions  and  two  different  forcing  functions. 

The  two  variable  coefficient  equations  are  solved  for  the 
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same  two  sets  of  boundary  conditions  and  for  one  forcing 
function.  The  problem  set  is  described  in  greater  detail 
in  Chapter  II. 

The  standard  finite  difference  method  is  used  only 
to  solve  the  linear  problems  whereas  the  finite  element 
methods  are  used  for  all  of  the  problems.  The  finite  dif¬ 
ference  formulations  for  each  problem  are  presented  in 
Appendix  B.  The  finite  element  formulations  for  each  prob¬ 
lem  are  described  in  Appendices  C  and  D. 

The  performance  criteria  chosen  are:  ease  of  formu¬ 
lation,  accuracy  of  the  solution,  computational  time 
required,  and  "efficiency"  of  the  numerical  method.  Effi¬ 
ciency  was  chosen  in  this  study  as  a  function  of  formation 
time  versus  accuracy  achieved;  where  accuracy  was  defined 
as  the  maximum  error  at  any  node  or  grid  point.  The  imple¬ 
mentation  of  the  methods  and  the  evaluation  of  their  perform¬ 
ance  was  done  with  the  use  of  computer  programs  written  by 
the  author  in  FORTRAN  EXTENDED.  Several  routines  were  used 
from  the  International  Mathematical  and  Statistical  Library 
(Ref  7)  for  solving  the  systems  of  equations  and  for  matrix 
manipulations.  The  programs  were  executed  on  the  Aero¬ 
nautical  Systems  Division  Computer  at  Wright-Patterson  Air 
Force  Base,  Ohio. 
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II.  Problem  Set  From  One-Dimensional 
Elliptic  Differential  Equations 

Second- Order  Linear — Axial  Rod 

The  axial  deflection  of  a  rod  under  the  influence 
of  an  axially  distributed  force  is  illustrated  in  Fig.  1. 
This  system  is  governed  by  the  differential  equation: 


_d_ 

dx 


|  A  (x)  E 


du(x)  | 
dx  j 


-F(x) 


(6) 


where 

E  =  Young's  modulus, 
u(x)  =  axial  deflection, 

F(x)  =  axially  distributed  force,  and 
A(x)  =  cross-sectional  area. 

For  a  tapered  rod  with  a  taper  ratio  a,  the  cross-sectional 
area  is  given  by 

A(x)  =  A  (1-aj)  ,  0<x<L  (7) 

where  A  is  the  area  at  the  root  of  the  rod.  If  the  taper 
ratio  a  is  zero,  the  area  is  constant  and  the  problem 
reduces  to  the  uniform  rod.  Two  possible  boundary  condi¬ 


tions  at  the  ends  are: 


(a)  Clamped-clamped,  half-sin  load,  ct=0 


(b)  Clamped-f ree,  uniform  load,  a=0.5 


Pig.  1.  Uniaxial  Deflection  of  a  Rod 
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u(x^)  =0  -  zero  displacement  (8) 


A(x)E 


du(x) 


dx 


=  0 


x=x, 


-  zero  axial  force  (9) 


where  x^  is  the  coordinate  at  the  particular  end.  The  six 
cases  of  the  axial  rod  studied  are  detailed  in  Table  I, 
where  the  uniform  load  is  F(x)  =  P  (constant);  and  the "half- 
sin"  load  is  P  (x)  =  P  sin  (-^)  . 


Fourth-Order  Linear — Beam  in  Bending 

The  governing  differential  equation  for  the  bending 
displacement  of  a  beam  subjected  to  a  transverse  load,  as 
illustrated  in  Fig.  2  is  fourth-order  and  linear.  The 
equation  is 


_dl 

dx2 


X  (x)  E 


d2w(x) 

dx2 


*  F  (x) 


(10) 


where  w(x)  is  the  transverse  displacement  and  I(x)  is  the 
cross-sectional  moment  of  inertia.  The  moment  of  inertia 
is  given  by: 

I(x)  =I(l-a-£-),  0<x<L  (11) 

where  I  is  the  moment  of  inertia  at  the  root  of  the  beam. 
The  boundary  conditions  studied  with  Eq  (10)  represent  the 
same  end  constraints  as  studied  with  the  rod.  The  boundary 
conditions  for  a  clamped  end  are: 
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(b)  Clamped- free,  uniform  load,  ct=0.5 


Fig.  2.  Transverse  Bending  Displacement  of  a  Beam 
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(12) 


w(x^)  =  0  -  zero  displacement 


dw 

dx 


=  0  -  zero  slope 


x=x. 


For  a  free  end  the  boundary  conditions  are 
d2w| 


I  (x)  E 


I  (x)E 


dx' 


d3w 

dx3 


=  0  -  zero  bending  moment 


x=x, 


(13) 


=  0  -  zero  shear  force 


x= 


*k 


As  with  the  axial  rod,  six  cases  are  studied.  In  all  cases 
the  boundary  conditions  at  the  left  end  are  given  by  Eq  (12) 
with  the  coordinate  xk  =  0.  The  details  are  given  in 
Table  II. 


Table  II 

Beam  Problems — Exact  Solutions 


Taper 

BC  Eqn 

Exact  Sol 

Description 

Ratio 

at 

Eqn  No 

Clamped-C lamped  Uniform 

(12) 

Beam,  Uniform  Load 

I 


Table  II — Continued 


Description 


Clamped-Clamped  Variable 
Beam,  Uniform  Load 


4LZ  \  PLx3 


Taper 

Ratio 


Exact  Sol 
Eqn  No 


w(x)  =  -gj-  J-  — -  +  —  (Cl-PL)  +  [  (2L-x)  In(l-^j-)  +  x]  (PIrCl  +  g4 


Cl  =  PL 


1 5  11  ,  n  ~  +  4  In  0.5 

|7  +  Tto0.5  (  #  C2  .  ^  24 - 3 - 


1  +  1.5  In  0.5 


1  +  1.5  In  0.5 


Taper 

Ratio 


Description 


Clamped-Free  Variable 
Beam,  Uniform  Load 


w(x)  =|j  L2[(2L-x)  In  (1--^)  +  x] 


Exact  Sol 
Eqn  No 
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Second-Order  Nonlinear 

A  one-dimensional,  steady,  second-order,  nonlinear 
analog  of  the  governing  differential  equation  for  transonic 
flow  over  an  airfoil  is  given  by 

(l-4#x)4)#xx  =  0  0<x<l  (20) 


where 


<)>  (x)  »  field  variable  -  potential, 

<j>,x  =  first  derivative  with  respect  to  x  of  the 
field  variable  <J>  ,  and 

4> ,  =  second  derivative  with  respect  to  x  of  the 

field  variable  <j> . 


The  boundary  conditions  are 
<Mo)  =  C1 
4>/x(°)  =  C2 


(21) 


and  either 


<J>  (1)  =  c4 


(22) 


or 

♦,x(l)  =  C3  (23) 

where  C^,  C 2,  C^,  and  are  constants. 


If  the  differential  equation  is  solved  with  the  two  boundary 
conditions  at  x=0  and  with  either  boundary  condition  at 
x=l,  the  exact  solution  may  be  discontinuous.  This  solu¬ 
tion  is  illustrated  in  Fig.  3  for  positive  constants 
through  C^. 
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<D(x) 

C4' 

■ 

_ ^ 

C1 

0 

4  V 

Fig.  3.  Exact  Solution  to  the 
General  Nonlinear  Problem 


The  coordinate  of  the  discontinuity  n  is  given  by 


C.+C-j-C . 
13  4 

C3"C2 


VC3 


The  solution  to  the  left  of  the  discontinuity  is 

4>l(x)  =  C2x+C1  ,  0<x<n  (25 

and  to  the  right  of  the  discontinuity  the  solution  is 

4>r(x)  =  C4(1-C3X)  ,  n-X-1  (26 

The  particular  problem  studied  in  this  work  has  constants 
C^=C4=0  and  C2=-C3=l,  thus  the  break  occurs  at  the  coordin 
ate  n=0.5.  This  solution  is  illustrated  in  Fig.  4. 

The  equations  of  the  solution  are: 


III.  Computer  System  and  Program  Information 

The  use  of  a  high  speed  digital  computer  is  a  must 
for  any  study  of  numerical  analysis  methods.  The  numerical 
calculations  for  the  linear  problems  studied  in  this 
thesis  were  performed  with  the  aid  of  the  ASD  computer  sys¬ 
tem.  This  system  has  two  Control  Data  Corporation  central 
processors  operating  in  parallel.  The  CDC  6613  and  CDC 
CYBER  74  processors  both  have  131000^  60-bit  words  of 
central  memory. 

The  computer  programs  written  by  the  author  for  this 
study  perform  all  but  three  of  the  manipulations  required 
to  formulate,  solve,  and  analyze  the  numerical  solutions. 
These  three  manipulations  are  performed  by  subroutines  from 
the  IMSL  code  in  operation  on  the  ASD  computer.  The  manipu¬ 
lations  and  associated  subroutines  are: 

1.  matrix  multiplication  —  VMULBB 

2.  linear  equation  solving  —  LEQT1B 

3.  vector  maximum  value  search  —  VABMXF 

All  matrices  are  stored  in  band  storage  mode  where  only 
elements  on  the  diagonals  are  stored.  A  flow  chart  of  the 
main  program  is  included  as  Fig  5. 
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IV.  Results  and  Discussion 

I 

The  results  of  this  study  are  analyzed  according 
to  the  types  of  differential  equations  solved.  Thus  the 
results  show  which  method  performs  best  for  each  class  of 
differential  equations. 

In  addition  to  the  performance  criteria  introduced 
in  Chapter  I,  several  other  criteria  are  also  included  in 
the  results,  for  they  show  a  more  complete  trend  than  just 
the  criteria  stated  in  Chapter  I.  The  entire  list  of 
criteria  is  given  in  Table  III.  When  comparing  errors 
the  roundoff  error  must  be  considered.  Because  of  round¬ 
off  error  in  the  solutions,  some  entries  in  Tables  A-I 
through  A-XII  are  nonzero,  but  many  orders  of  magnitude 
smaller  than  similar  entries.  The  growth  in  these  nonzero 
terms  is  of  great  importance  for  they  indicate  an  unstable 
algorithm  when  the  growth  rate  is  exceptionally  large.  This 
phenomenon  is  discussed  later.  In  all  cases,  the  error  com¬ 
pared  is  a  relative  error  eR  defined  at  a  point  as 


e 


R 


Hx^T 


where 

$  *=  approximate  solution  at  a  point  in  the  domain, 

$  ■  exact  solution  at  the  same  point  in  the  domain, 
and 
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Table  III 


Performance  Criteria  Used  for  Comparing  the 
Numerical  Methods  for  Linear  Problems 


Table  Entry 

Description  of  Criteria 

Max  Nodal 
Error 

— 

Maximum  relative  error  at  any  node  or 
grid  point 

Max  Point 
Error 

Maximum  relative  error  at  any  one  of 

97  equally  spaced  points  in  the 
domain — all  nodes  are  included  in 
this  sample 

Best  Nodal 
Error 

— 

Best  Maximum  Nodal  Error  for  each 
method  and  each  problem 

Best  Point 
Error 

— 

Best  Maximum  Point  Error  for  each 
method  and  each  problem 

Execution 

Time 

— 

Decimal  seconds  required  to  solve 
the  global  equations 

Formation 

Time 

Decimal  seconds  required  to  compute 
and  assemble  the  global  equations; 
and  implement  the  boundary  conditions 

CF 

— 

Convergence  factor  for  a  particular 
halving  of  the  stepsize 

CF  e (h) 

CF  "  e (h/2) 

— the  order  of  convergence  is  equal 
to  v^2f.  If  CF<1.0  the  algorithm 
diverges 

CFN 

— 

Max  Nodal  Error  convergence  factor 

CF  for  smallest  stepsize  studied 

CFP 

— 

Max  Point  Error  Convergence  factor 

CF  for  smallest  stepsize  studied 

Efficiency 

Point 

Maximum  efficiency  for  each  method 
and  each  problem  where  efficiency  is 
the  log  f [Max  Point  Error  (Formation 
Time) 5] “1} 

Number  of 
Elements 

— — * 

Number  of  Elements  corresponding  to 
maximum  efficiency 

(x) .  =  exact  solution  at  Xi=0.5  for  clamped- 
clamped  boundary  conditions  and  X.=1.0 
for  clamped-free  conditions. 

These  points  give  an  exact  solution  that  is  either  the 
maximum  or  very  near  the  maximum  displacement. 

Second-Order  Linear — Axial  Rod 

The  results  presented  in  Tables  IV  through  IX  indi¬ 
cate  a  general  trend.  That  is,  the  higher-order  finite 
element  methods,  collocation  and  least  squares,  are  more 
accurate  throughout  the  domain  than  are  finite  differ¬ 
ences  and  Galerkin's  method.  This  is  evidenced  by  small 
pointwise  errors  for  the  higher-order  methods.  Least 
squares  is  by  far  the  most  accurate  and  most  efficient 
method  studied  for  this  second-order  problem. 

Finite  Differences  and  Galerkin's.  These  two 
methods  both  use  a  linear  approximation  for  the  field  vari¬ 
able.  This  order  of  approximation  is  not  as  accurate  as 
the  cubics  even  though  the  solution  at  the  nodes  matches 
the  exact  solution  identically.  For  Problems  1  and  2, 
the  values  recorded  in  Table  IV  and  V  for  these  two  methods 
correspond  to  a  solution  with  thirty-nodes.  The  Best  Point 
Error  indicates  the  maximum  error  at  one  of  97  equally 
spaced  points  in  the  domain.  This  quantity  reveals  an 
accuracy  of  only  two  digits  for  these  methods.  Both 
methods  proved  to  be  accurate  to  the  order  h2  as  expected, 
for  all  the  second-order  problems  except  for  Problem  1 
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Table  IV 


Results — Problem  3?  Variable  Clamped-C lamped  Rod,  Uniform  Load,  a=0.5 


where  the  nodal  error  increased  at  a  rate  proportional  to 

h-/576  =  h“2.4  _  According  to  Fix  (Ref  2:205-215), 

the  round-off  error  is  bound  by  the  condition  number  of  the 

Toeplitz  matrix.  This  upper  bound  is  h~  where  m  is  the 

order  of  the  differential  equation.  In  this  case,  the  round- 

-4 

off  error  should  be  less  than  Ch  ;  C=constant.  Therefore, 
the  errors  present  in  these  problems  are  acceptable  con¬ 
sidering  that  the  error  present  is  in  the  thirteenth  digit 
on  a  machine  with  fourteen  significant  digits. 

Collocation  and  Least  Squares.  These  two  methods 
give  excellent  accuracy  for  Problems  1  and  2  but  are 
definitely  unstable.  This  is  due  to  the  fact  that  the  cubic 
approximation  with  two  boundary  conditions  will  solve  the 
governing  differential  equation  for  Problems  1  and  2  exactly 
with  only  one  element.  By  increasing  the  number  of  ele¬ 
ments  one  simply  increases  the  probability  that  the  dis¬ 
cretization  error  and/or  the  roundoff  error  will  be 
increased.  On  the  other  hand,  for  the  four  remaining  rod 
problems  the  exact  solutions  are  either  functions  of  natural 
logarithms  or  trigonometric  functions.  These  solutions  can 
not  be  modeled  exactly  with  a  third-order  polynomial.  For 
these  four  problems  the  collocation  method  is  only  of  order 
h2  where  as  the  least  squares  method  approaches  an  order 
of  accuracy  h*.  It  is  also  interesting  to  note  that  the 
order  of  the  approximation  polynomial  does  not  noticeably 
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affect  the  round-off  error.  This  can  be  seen  from  the 
results  of  Problems  1  and  2,  where  for  the  collocation 
method  the  error  grows  at  approximately  the  same  rate  as 
for  the  finite  difference  and  Galerkin  methods.  This 
observation  is  also  made  by  Fix  (Ref  2  :215). 

Fourth -Order  Linear- -Beam  in  Bending 

The  results  for  this  system  are  given  in  Tables  X 
through  XV.  The  tables  show  that  for  the  simple  cases  of 
Problems  7  and  8  the  methods  using  a  septic  approximation 
are  exceptionally  accurate  and  very  efficient.  But  as  the 
problem  becomes  more  complex  the  septic  shape  functions  do 
not  guarantee  a  more  accurate  solution  than  is  obtained  by 
the  Galerkin' s  method  with  cubic  shape  functions.  It  is 
also  apparent  that  the  order  of  convergence  is  very 
strongly  dependent  on  the  particular  problem  being  solved. 
This  fact  is  substantiated  by  examining  the  convergence 
factors  CF  for  the  Galerkin  method  in  Tables  X  through  XV. 
For  the  tapered  beam  the  algorithm  is  only  of  order  h2 
whereas  for  the  uniform  beam  with  a  half- sin  load  the 
method  is  of  order  h4 .  In  Problems  7  and  8  the  change  from 
0(h4)  to  0(h2)  is  brought  about  by  simply  freeing  the  right 
end  of  the  beam.  As  with  the  Axial  Rod,  the  higher-order 
methods  are  unstable  for  the  cases  where  the  exact  solution 
is  of  less  order  than  the  approximation  function.  The 
explanation  of  this  is  the  same  as  for  the  axial  rod.  For 
the  simple  cases  the  collocation  or  least  squares  methods 


Results — Problem  7;  Uniform  Clamped-Clamped  Beam,  Uniform  Load 


001  . 906E-11  .001  . 906E-11  15.08 


Table  XIV 


048  . 601E-14  15.7  .535E- 


give  good  results  with  only  one  element.  The  collocation 
method  is  the  most  accurate  and  efficient  method  for  all  of 
the  problems  except  the  half-sin  loaded  beams.  For  these 
last  two  problems  the  Galerkin  method  was  more  efficient 
due  to  the  more  simple  calculations  required  to  form  the 
global  equations.  The  collocation  me '-hod  exhibited  very 
poor  convergence  for  the  variable  beams. 

Nonlinear  Problem 

The  nonlinear  problem  was  solved  without  the  aid  of 
the  digital  computer  due  to  a  lack  of  time.  The  two  approxi¬ 
mate  methods  used  were  the  Galerkin  method  and  the  colloca¬ 
tion  method.  Two  iterative  formulations  of  the  original 
differential  equation  were  used  as  detailed  in  Appendix  F; 
each  proved  to  behave  quite  differently.  It  was  expected 
that  formulation  A,  which  does  not  have  an  iterative  forcing 
term,  would  give  the  fastest  convergence.  This  was,  in  fact 
the  case;  formulation  A  gave  the  exact  solution  after  one 
iteration.  Formulation  B  converged  to  the  correct  solution 
but  took  several  iterations.  The  numerical  results  for  this 
problem  by  the  Galerkin  method  are  given  in  Table  XVI, 
and  for  the  collocation  method  in  Table  XVII.  The  process 
of  solving  the  equations  using  a  hand  calculator  forces  one 
to  place  more  emphasis  on  the  "efficiency"  of  the  method. 
Efficiency  in  this  context  could  be  defined  as  the  number 
of  keystrokes  required  to  obtain  the  solution  for  each 
iteration.  By  this  measure  the  Galerkin  technique  is 
definitely  easier  to  work  with. 
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Table  XVI 


XVII 


The  nonlinear  problem  studied  required  some  knowl¬ 
edge  of  the  physical  process  that  Eq  (F-l)  is  used  as  an 
analog  for.  The  boundary  conditions  chosen  are  mathe¬ 
matically  convenient  but  not  physically  desirable.  The 
condition  <j>(l)=0  implies  that  the  function  itself  goes  to 
zero  at  x=l.  A  better  choice  for  the  boundary  condition  at 
x=l  would  have  been  <J>X(1)=C^.  This  implies  that  the  func¬ 
tion  is  nonzero  at  the  boundary.  Another  aspect  of  the 
problem  which  must  be  understood  is  the  condition  that  a 
jump  in  the  first  derivative  of  <p  must  occur  in  the  flow 
field.  The  function  is  continuous  but  the  derivative 
changes  discontinuously /  which  models  the  idealized  behavior 
caused  by  a  shock  in  the  flow  field.  This  behavior  requires 
an  approximation  to  the  field  variable  which  has  first 
derivatives  as  nodal  parameters.  Such  an  approximation  dis¬ 
qualifies  linear  shape  functions  even  though  the  exact  solu¬ 
tion  can  be  represented  by  two  straight  lines. 
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V.  Conclusions 


One  of  the  objectives  of  the  study  was  to  determine 
if  comparisons  made  by  Houstis  (Ref  5)  based  on  a  second- 
order,  two  dimensional,  constant  coefficient  problem  set 
would  hold  for  a  set  of  higher-order  variable  coefficient 
problems  or  nonlinear  problems.  In  searching  for  an  answer 
to  this  question  several  other  conclusions  were  also 
obtained.  The  conclusions  of  this  study  are: 

1.  For  the  second-order  linear  problem  studied, 
the  least  squares  method,  with  Hermite  cubic  polynomials, 
is  more  accurate  and  also  more  efficient  than  the  colloca¬ 
tion,  Galerkin,  or  finite  difference  methods  studied.  The 
least  squares  method  outperformed  the  other  four  methods  for 
all  six  cases  of  the  axial  rod  studied.  The  most  heavily 
weighted  performance  criteria  were  efficiency  and  accuracy, 
in  that  order. 

2.  For  the  fourth-order  linear  problem  studied, 
the  least  squares  method  with  Hermite  septic  polynomials 
was  only  used  for  the  constant  coefficient  case  with  a  con¬ 
stant  forcing  term.  In  this  case  least  squares  performed 
equally  as  well  as  collocation,  with  the  same  interpolation 
functions;  and  significantly  better  than  finite  differences 
and  Galerkin  with  lower-order  interpolation  functions. 

For  the  variable  coefficient  case  collocation  was  slightly 
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less  efficient  than  the  other  two  methods.  In  general, 
collocation  was  the  method  of  choice. 

3.  For  the  nonlinear  problem  the  Galerkin  and 
collocation  methods  both  perform  similarly.  The  number  of 
calculations  required  for  each  iteration  was  much  greater 
for  the  collocation  method  than  for  the  Galerkin  method. 

The  second  conclusion  obtained  by  Houstis  varies 
significantly  from  those  of  this  study.  This  difference 
can  possibly  be  accounted  for  by  the  fact  that  Houstis  used 
a  numerical  integration  scheme  for  the  Galerkin  and  least 
squares  methods.  Since  the  methods  were  formulated  with 
polynomials,  the  integration  is  trivial  and  numerical  inte¬ 
gration  is  not  called  for. 
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Numerical  and  Graphical  Data 

This  appendix  contains  all  of  the  numerical  results 
for  the  linear  problems.  These  results  are  given  in 
Tables  A-I  through  A-XII  for  which  the  headings  are  defined 
in  Chapter  IV.  Each  table  presents  data  for  one  specific 
problem  as  solved  using  the  various  numerical  methods. 

The  graphical  data  shows  two  aspects  of  the  study. 
The  first  six  graphs  are  representative  of  the  way  in  which 
the  methods  approximate  the  solutions  to  the  various  prob¬ 
lems.  The  next  twenty-four  graphs  are  plots  of  the  relative 
error  at  97  points  in  the  domain.  This  number  of  points 
was  chosen  because  it  insures  that  the  plotted  solution  will 
not  "skip  over"  any  nodes  for  the  values  of  E  studied. 


Table  A- I 


Uniform  Clamped-Free  Rod#  Uniform  Load 


Table 
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Table  A-III 


Table  A-III — Continued 


Table  A-IV 


Table  A-IV — Continued 
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Table  A-V — Continued 


Table  A- VI 
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Table  A- VI I — Continued 


Table  A-VIII 


Table  A-VIII — Continued 


Table  A- IX 
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Table  A-IX — Continued 


Variable  Clamped-Free  Beam,  Uniform  Load,  Taper 


Table  A-X — Continued 


A-XI 


Table  A-XI — Continued 
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ble  A-XII — Continued 
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RELATIVE  ERROR 
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Appendix  B 

Finite  Difference  Expressions  and  Sample  Problems 

The  finite  difference  method  uses  difference  quo¬ 
tients  to  approximate  the  derivatives  in  a  differential 
equation.  These  quotients  are  formed  at  N+l  (a  finite 
number)  grid  points  in  the  problem  domain.  For  the  one¬ 
dimensional  problems  in  this  study  the  domain  is  discretized 
as  shown  in  Fig  B-l. 


Fig  B-l.  Discretization  of  One-Dimensional  Domain 
for  Finite  Difference  Method 

Near  the  boundaries  the  difference  quotients  may  require 
that  a  "false"  node  exist  outside  the  domain.  This  situa¬ 
tion  can  easily  be  handled  if  symmetry  about  the  boundary 
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r 


is  assumed  (Ref  8:143).  This  situation  is  illustrated  in 
Pig  B-2. 


Fig  B-2.  Symmetry  of  Approximate  Displacement 
About  a  Clamped  End 

Axial  Rod 

The  governing  differential  equation  for  the  axial 

rod  is 

A{a(x)  E  $i}  =  -Fix)  (B-l) 

which  after  differentiation  becomes 

A(x)E  —  +  E  =  -F(x)  (B-2) 

du2  dX  035 

The  two  second-order  central-difference  quotients  needed 
are 
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cMvi)  '  7  (yj-i-2yj+yj+i>  +  <-AhVv> 


(B-3) 


D  (yi)  -  ^  (-Yj.!  +yj+1)  +  (B 

Substituting  Eqs  (B-3)  and  (B-4)  into  Eq  (B-2)  gives  the 
result 


[— A(x j__1)  -4A(Xj)  +  A(Xj+1)]  +  8A(Xj)Uj 

4F(x.)h2 

+  [A(Xj_1)-4A(Xj)  -A(Xj+1)  )Uj+1® - - 

j=0,l,...,N  (B-5) 


For  a  uniform  rod  Eq  (B-5)  reduces  to 


-  uj-i  +  2uj  “  uj+l 


F  (x..)h2 
EA 


j=0,l,...,N  (B-6) 


where  F(x^)  »  P  for  the  uniform  load  case  and  F(x..)  = 
irx . 

P  sin  (-s-J-)  for  the  half-sin  load  distribution,  and 

Jj 

x-i 

A(x.)  *  (l-o -r^  ).  In  this  study  only  two  boundary  condi- 
3  ^ 

tions  are  studied,  they  are: 


u(xk)  a  o  -  zero  displacement  (B-7) 
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(B-8) 


A(VE  i 


(xk)  = 


0  -  zero  axial  force 


To  account  for  the  boundary  conditions,  the  equa¬ 
tions  associated  with  nodes  near  the  end  of  the  rod  are 
modified  as  follows: 

1.  For  a  clamped  end  at  xk=0,  then  uQ=0  in  row  1 
and  in  row  2 

4F (x, )  h2 

+  fA(xQ)-4A(x1)-A(x2)  ]ui2  - - -  (B-9) 


2.  For  a  clamped  end  at  Xk=L,  then  ^=0  in  row 
N+l  and  in  row  N 


8A(xN-1)uN-1  +  [A(xN-1)“4A(xN-2)_A(xn-3)]uN-2  “ 


4F(xN-l)h2 

E 


(B-10) 


3.  For  a  free  end  at  x^=L,  then  in  row  N+l 

4F(xM)h2 

8a<V“n  -  8a<*n-i)'Vi - -r—  «B-n» 


Sample  Problem  B-l .  A  clamped-free  uniform  rod  is 
loaded  with  a  half-sin  distribution  of  uniaxial  tension. 

The  rod  is  discretized  with  five  grid  points.  After  the 
boundary  conditions  are  applied,  the  system  equations  become 
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The  approximate  solution  compared  with  the  exact  solution 


approximate 


exact 


.15088 


.25758 


.30178 

.30178 


.15122 


.26048 


.31038 

.31831 


Beam  in  Bendinc 


The  governing  differential  equation  for  the  beam 


in  bending  is 


—  1 1  (x)  E  —  1  =  F  (x) 

r»v2  '  Av2  ‘ 


(B-12 ) 


After  differentiation  it  becomes 


X(x)  d4w  +  2  dl (x)  d3w  +  d2I (x)  d2w  _ 
dx4  dx  dx3  dx2  dx2 


(B-13) 


The  fourth-order  central-difference  quotients  needed  are 

D<yJ>  -  T5h(yj-2-9yj-l+8yj+ryj+2l  +  <T0  h'yjV' 
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(B-15) 


Dty  =^?'-yj-2+16y3-l-30y:+16Y^'y3«)  +  (^h'yjVI) 


DS(y3’  '^r('y3-2+2y3-i'2y3^j«)  +  (-7h!yjV)  <B'16> 

°* <y3>  -  ffr  0V2-‘yH«y3%1%2l  +  <-?  h'y™) 


If  the  moment  of  inertia  I  (x)  is  assumed  to  be  a  constant 
or  a  linear  function  then  Eq  (B-14)  simplifies  to 


Kx) 


d4w 

dx4 


+  2 


dl (x)  d3 


dx 


dx 


w  _  F  (x) 
3  E 


Substituting  the  central-difference  expressions  for  the 
derivatives  gives 


1<xj_l)wj_2  [! (Xj_^)+I (Xj) ]w^_^  + 


[I  (x  +41  (x^)  +1  (x  j+1)  ]Wj 


-2[I(Xj)+I(Xj+1)]Wj+1  + 


I*xj+2*wi+2 


=  E(x)h4 
E 


j=0,l,. . .  ,N 


where  F(x.)  =  P  (constant)  for  the  uniform  load  distribu- 
3  ttx  . 

tion  and  F(x.)  =  P  sin  (—=-*)  for  the  half- sin  load  dis- 
D  b 

tribution. 

The  boundary  conditions  modify  the  set  of  system 
equations  as  follows: 

1.  For  a  clamped  end  at  x^=0,  then  Wq=0  in  row  1, 

and 
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' 


r  I 


and 


[21 (x0)+4I(x1)+I(x2) ]w1-2 [I (x1)+: (x2) ]w2 


F(Xl)h" 

+  I(x2)w3  =  - g -  in  row  2,  and 


-2  [I  <x0)+I  (Xj^)  ]w1+[I  (x2>+41  (x2)+I  (x3)  ]w2 


F(x2)h*' 

•2[I(x2)+I  (x3)  ]w3+I  (x3)w4  =  - g -  in  row  3 


2.  For  a  clamped  end  at  -x^L,  then  wN=0  in  row  N+l, 


I<XN-2)wN-3-2 [I (xN-23+I (N-13 lwN-2 


F(xN_  )h- 

+  [2I(xN)+4I(xN_1)+I(xN^2)]wN_1= - Y -  in  row  N, 

and 

I(xN-3,WN-4-2fI^XN-2)+I (xN-33  3wN-3 
+  ^  ^^-1^  +41  ^-2 3  +I  ^^-3 3  3wN-2~2  ^  ^XN-23 

P(aWh* 

+  1(xN-l33wN-l  =  - E -  in  rOW  N_1 

3.  For  a  free  end  at  x^L,  then 

FfXjjJh4 

21  ,XN-l)"N-2-41  'Vi'Vi*21  ( VWN  =  - E - 

in  row  N+l,  and 
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1 (xN-2)wN-3~2[I (xN-2)+(xN-1) JwN-2 

+  II(xN-2)+4I(xN-l)lwN-l~  2I(xN-1)wN  in  row  N 

Sample  Problem  B-2 .  A  clamped-free  uniform  beam 
is  subjected  to  a  half -sin  transverse  load  distribution. 
The  beam  is  discretized  with  five  grid  points.  After  the 
boundary  conditions  are  applied,  the  system  equations 
become 


—  — 

/  \ 

f  \ 

7-4100 

w2 

0.002762 

-4  6-4  1  0 

w, 

0.003906 

< 

3 

►  =  i 

> 

1-3  6-4  1 

w4 

0.002762 

002-42 

w5 

^  > 

0.0 

The  approximate  solution  compared  with  the  exact  solution 
is 

approximate  exact 


> 

f 

f  > 

r 

w2 

0.00943 

w2 

0.00831 

i  W 3 

w4 

►  - 

PLl  < 

IE  1 

0.02829 

0.04991 

»  < 

W3 

W4 

IE  ' 

0.02729 

0.05021 

w,- 

0.07154 

W- 

0.07385 

D 

\  > 

; 

k  5  > 

J 
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Appendix  C 

Derivation  of  Finite  Element  Equations  for  the 
Axial  Rod  with  Sample  Problems 


The  governing  differential  equation  for  the  axial 
rod  under  the  influence  of  an  uniaxially  distributed  force 
is  given  as 


dx 


1 A  (x)  E 


du(x) 

dx 


=  -F  (x) 


(Ol) 


Discretizing  the  domain  into  E  elements  with  E+l  nodes, 
the  solution  can  be  approximated  as 


u (x)  =  N j (x>Uj  j=l,2,...,M 


where  M  =  number  of  nodal  parameters,  and  the  repeated 
index  implies  summation. 

The  differential  equation  will  not  be  satisfied 
exactly.  The  residual  error  e  is  thus  defined  as 


dN  . 


e(x)  =  £  jA(x)E  -jJ-  u.  (+  F (x)  *  0 

j-1,2, . ..  ,M 


(C-2) 


The  method  of  weighted  residuals  (Ref  12)  forces  the  residual 
error  to  be  zero,  in  some  weighted  average  sense,  over  the 
domain  fi.  For  each  method  this  is  achieved  by  performing 
the  following  integrations. 
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Galerkin: 


I 


e (x) 6 (x-x^) dx  = 


e (x)N^ (x) dx  =  0 


0 


Least  Squares: 


I 


9e  (x) 
3ui 


dx  =  0 


(C-3) 


(C— 4 ) 


(C-5) 


Collocation  Method 

The  evaluation  of  the  integral  for  the  collocation 
method  is  carried  out  by  first  discretizing  the  domain  into 
E  elements,  each  of  length  h.  The  integral  is  thus: 

E  /.eh 

.  I  =  yi  I  e (x) 6 (x-x^) dx  =  0  (C— 6 ) 

e=l  •'(e-l)h 

This  integral  can  be  simplified  by  using  local  coordinates 
x,  where 

x  =  (e-1) h+x ,  0<x£h  (C-7) 

0<_x£L 


,1 


3 


The  E  elemental  integrals  with  this  transformation  become 
very  similar;  that  is: 


-I' 


e.(x)  6  (x-x^) dx 


=  0 


(C-8) 


Substituting  Eq  (C-2)  into  Eq  (C-8)  along  with  the  coordiu 
ate  transformation  gives 


I 


e 


6 (x-x^) dx  =  0 


3*1,2, • • « ,m 
i“l , 2 , • • • , m 


(C-9) 


where  N .  (x)  are  the  local  shape  functions  (tabulated  in 
3 

Appendix  D)  and  m  is  the  number  of  nodal  parameters  per 
element . 

Integration  of  Eq  (C-9)  gives 


"e  * 


—  Ia(x) 

dx  L 


dN.(x)  -I 

E  — •*- -  u.  +F(x) 

dx 


*  0 


x*x. 


3*1 , 2 , • » • ,m 


(C-10) 


Suppose  the  area  A(x)  is  assumed  to  vary  linearly  in  the 
domain,  then  it  can  be  expressed  as 


A (x)  -  A  |l  -  j  t(e-l) 


(C-ll) 


104 


where 


A  =  area  of  element  at  left  end  (x=0) ,  and 
a  =  taper  ratio  (0£a<1.0). 


Substituting  Eq  (C-ll)  into  Eq  (C-10)  and  differentiating 
yields 


d2N . 


a  -  )  u 

?•[  (e-l)h+x]  ?AE  - 1 


u . 

dxa  3 


Y  AE 

1j 


dN  . 


un  .  \ 

— 1  u.+F(x)> 
dx  3  }\ 


=  0 


x=x. 


i  1 , 2  ,  •  •  •  ,m 
j =1 ,2 ,  •  • . ,m 

e=l,2, . . .  ,E  (C-12) 

d2N . 

where  N.  must  be  such  that  - ^  is  non-zero.  In  this  thesis 

3  d£2 

Nj  is  chosen  to  be  the  cubic  shape  functions  associated  with 
the  C1  line  element  (see  Table  E-I) .  For  this  class  of 
functions  there  are  four  nodal  parameters  per  element 
(m  equals  four) .  With  cubic  shape  functions  (Eq  (C-12)  can 
be  written  for  element  e  as 


Uj  +  F. 


0 


i=l, . . . , 4 
j=l # • • • t 4 
e  ~  1  , « .  •  ,  E 


(C-13) 


where 

K^j  =  elemental  stiffness  matrix, 

Uj  =  vector  of  element  nodal  parameters,  and 
F^  =  vector  of  element  equivalent  forces. 


In  the  local  coordinate  system,  the  elemental  stiff¬ 


ness  matrix  has  elements  of  the  form 


Ki2=f  I'-lMIt-414 


|je-l)  J  }[“2+6 


(C-14) 


where  i=l,2,3,4  and  the  x^  have  been  chosen  as  the  equally 
spaced  points  jL  =  [0,  h]  . 

The  choice  of  equally  spaced  points  over  Gaussian  points 
was  due  to  the  conclusion  of  Houstis  (Ref  5:327)  "that 
equally  spaced  points  give  slightly  better  accuracy  for 
rectangular  domains."  Theoretically  the  Gauss  points  give 
better  results  according  to  Prenter  (Ref  10:304-314).  He 

showed  that  "collocation  at  Gaussian  points  using  a  basis 

4 

of  piecewise  cubic  Hermite  polynomials  gives  an  order  h 
algorithm  for  approximating  the  unique  solution  x(t)  to  the 
equation 


-x"(t)  +  o  (t )  x  ( t)  =  f(t) 


a<t<b 


For  a  uniform  rod  the  taper  ratio  a  equals  zero. 
This  greatly  simplifies  the  stiffness  matrix. 

Thus: 


6 

-4h 

6 

-2h 

[K]  =  — 

-2 

-2h 

2 

0 

h 

2 

0 

-2 

2h 

6 

2h 

-6 

4h 

the  local  coordinate  , 

system, 

the 

local 

(C-15) 


are  related  to  the  global  unknowns  as  follows: 


local  global 


'  u 

(0) ' 

'  ue-l  ' 

u ( (e-1) h) ' 

U2 

< 

’  =  < 

du 

dx 

(0) 

►  =  < 

ue 

►  =  < 

|j((e-l)h) 

U3 

u 

(h) 

ue+l 

u(eh) 

'u4- 

du 

dx 

(h) 

,  ue+2  ' 

^(eh)  - 

(016) 


The  local  and  global  force  vectors  are  related  in  the  same 
manner. 

local  global 


F1 

F  (0) 

Fe-1 

F((e-l)h) 

F2 

F  (§) 

Fe 

F ( (e-1) h+h/3) 

►  =  i 

►  =  < 

y 

F3 

F  (f . 

Fe+1 

F( (e-l)h+2h/3) 

» F4  * 

F  (h) 

*  Fe+2 , 

F  (eh) 
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where  F(x)  =  F(x)  =  P  (constant)  for  the  case  of  uniform 
load;  and  F(x)  *  P  sin  <^~)  for  the  case  of  a  half-sin  load 
distribution. 

When  elements  are  connected  together  the  shape 
functions  chosen  guarantee  continuity  of  the  function  and 
its  first  derivative  in  the  interval  (0,L) .  Therefore, 
when  forming  the  global  equations,  the  local  equations  are 
added  so  that  their  local  nodal  parameters  coincide  with  the 
appropriate  global  parameters  as  defined  in  Eq  (C-16) . 

Sample  Problem  C-l .  A  clamped-free  uniform  rod 
with  a  half-sin  load  distribution  is  discretized  with  two 
line  elements.  The  nodal  parameters  are  as  illustrated 
in  Fig  C-l. 


Fig  C-l.  Global  Illustration  of  a  Clamped-Free  Uniform 
Rod  Approximated  by  Two  Cl  Line-Elements  with 
a  Half-Sin  Load  Distribution 
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The  unreduced  global  equations  are 
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-2L 

6 

-L 
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f  -s 

U1 

sin(0) 

+  0 

-2 

-L 

2 

0 

0 

0 

U2 
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U3 

►  =  P*< 

sin(-j-) 

+  Sin(-J-) 
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6 

L 

-8 

L 

2 

0 

u4 

sin(^-) 

+  sin(yO 

0 

0 

2 

0 

-2 

L 

U5 

0 

+  sin(^p) 

0 

0 

6 

L 

-6 

2L 

^  U6> 

0 

+  sin(ir) 

These  equations  can  be  reduced  by  using  the  boundary  condi¬ 
tions.  At  the  left  end,  the  displacement  u^  is  zero;  and 
at  the  right  end,  the  axial  force  AE  ug  is  zero.  Therefore 
rows  one  and  six  and  the  corresponding  columns  can  be 
eliminated  from  the  system  of  equations.  The  result  is 


-L  2  0  0 

f 

u? 

1 

,  ' 

.5 

0  4  -L  6 

1.866 

< 

J 

►  =  P  « 

L  -8  L  2 

U4 
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0  2  0  -2 

uc 

V  5 

>  *5  . 

The  approximate  solution  compared  with  the  exact  is 
approximate  exact 


u2 

.64431 

’  U2 

r  .6366' 

u3 

PL2 
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u3 
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u4 

’  *  AE  4 
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»  Up  4 
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Galerkin  Method 


Galerkin's  Method  is  expressed  by  Eq  (C-4) .  That 

is: 

I-f  {ak  [A(X,E  +  p(=‘>}  Vx  =  o 
J0 

i=l f  2  f • • • t M 


Integration  by  parts  of  the  first  term 


dN.  dN . 


A(*>E  I?  ~Sx  V*  +  A<*>E  S  Ni 


{ 


L  .1> 
+ 

0 


F(x)Nidx  =  0 


i— 1 , 2 , . . . ,M 

j=l r 2 , . . . ,M  (C-18) 


du 


The  boundary  term  A(x)E  ^  tT 


vanishes  for  both  cases  of 


|0 

boundary  conditions  studied  in  this  thesis.  Therefore 
Eq  (C-18)  becomes 


I 


A  (x)  E 


dNi 

dx 


dx 


u . 
3 


F(x)N^dx  =  0 


i=l r 2 / 
j=l,2, 


,M 

|M 


Since  the  domain  is  discretized  into  E  elements  the  inte¬ 
gral  I  can  be  written  as 
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For  these  integrals  to  exist  the  shape  functions  must 
be  at  least  C°  continuous.  They  are  chosen  in  this  study 
to  be  the  linear  shape  functions  associated  with  the  C° 


>  line  element  described  in  Appendix  E.  For  the  Galerkin 

r 

method  Eq  (C-ll)  is  more  conveniently  expressed  as 
I  A(x)  «  (1  -  jj)  +  (|)Ar  (C-20) 

I 

where 

\  -  »[i  - 

*E  '  A[X  -  ¥\ 

,  I 

Equation  (C-7)  can  be  used  to  transform  Eq  (£-19)  to  a 
t  local  coordinate  system.  With  the  substitution  of 

Eq  (C-20)  the  element  integral  Ie  becomes 


e=l,2,...,E  (C-21) 
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where  F(x)  =  P  (constant)  for  the  uniform  load  distribution; 
and  F(x)  =  P  sin  (£[e-l) h+  x] )  for  the  half-sin  load  dis¬ 
tribution.  For  the  half-sin  load  distribution  Eq  (C-21) 
becomes 


I 


e 


<W  = 

2h 


(ire 
sm  — 

E 


+ 


PL 

IT 


-cos 


TT  (e-1) 
E 


cos 


ire  l 

E  ' 


0 

0 


e=l, 2 , ,E  (C-22) 


For  the  uniform  load  case  and  constant  area  A,  Eq  (C-21) 
becomes 


Je  - 


AE 


'  1 

-i 

ue 

Ph 

1 

= 

0 

r1 

i 

ue+l 

2 

1 

0 

e=l, 2 , . . . ,E  (C-23) 


These  equations  are  identical  to  the  equations  for  the 
standard  finite  difference  method  detailed  in  Appendix  B. 

Sample  Problem  C-2.  A  constant  area  clamped-free 
axial  rod  is  loaded  with  a  half-sin  distribution  of  uni¬ 
axial  force.  This  problem  is  approximated  with  two  C°  line 
elements  as  shown  in  Fig  C-2. 
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which  matches  the  exact  solution. 


Least  Squares  Method 

The  least  squares  technique  of  forming  finite  ele¬ 
ment  equations  relies  on  the  evaluation  of  Eq  (C-5) .  Sub¬ 
stituting  Eq  (C-2)  into  Eq  (C-5)  yields 

1  uj]  +  F,x>j  ah  [R(X)E  - 0 

i=l , 2  , « . • , M 

j— 1,2 , , . . ,M  (C-24) 


Multiplying  and  differentiating  terms  in  the  integrand 
gives 


I  =  E 


'/.  I*' 


d2N. 


d2N_ 


(x) 


+  2 


dx2  dx"1 


dA(x) 

dx 


dN. 
_ l 

dx 


dN 


-a  J  dx  u . 
dx  j 


F(X) 


d2N. 

A(x) - 1 

dx 


+ 


dA  (x) 
dx 


0 


i=l ,2,...,M 
j— 1 ,2,...,M 


(C-25) 


This  integral  can  be  written  as  the  sum  of  the  ele¬ 
mental  integrals  in  a  local  coordinate  system.  The  trans¬ 
formation  from  global  to  local  is  given  by  Eq  (C-7) .  The 
global  shape  functions  could  also  be  transformed,  but  the 
linear  shape  functions  of  Table  E-l  have  all  the  required 
properties  of  the  local  shape  functions  desired.  To 
simplify  the  albegra  the  derivatives  of  the  shape  functions 
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are  written  as 


^  ■ h  [v*i  (I)  ♦  x  (!)  ] 


i=l,2, ... ,4 


d2N. 

- -  =  A. 

dx2  1 


1=1, 2,.  ..,4 


where 

ai  =  [0'  P  °'  01 

A.  =  [-6,  -4h,  6,  -2h] 

1  h2 

B.  =  —  [12,  6h,  -12,  6h] 

1  h2 

The  final  result  after  integrating  the  polynomial  terms 
with  F(x)  =  P  (constant)  can  be  written  for  element  e  in 
the  form 


I  =  K.  -u  . 
e  13  3 


+  F±  =  0 


l—l ,2,...  ,4 
3=1 r 2 , . • .  ,4 
e=l,2 , . . . , E  (C— 26) 


where 

Kij  '  E 

+  '•iV'iV'WV' 

+  («iYBi“j,|-KsR,5Vi 
+  AiAjt?5  (20ai.!  '  40  alV80V)1 
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I 


+  (A.Bj+BiAj)[i  (5AL!-20ALAR+45AB!|] 

+BiBj<^  ,2V-9alV27V>'} 


i=l,2,3,4 

j«l,2,3,4 


Fi  =  p  'WV  +  (Ai  +  I  Bi»V 


*  (e-1) 


du  it  (e-1) 
dx  - 


E 

du  ire 
dx  - 


For  a  uniform  rod  the  elemental  stiffness  matrix  is 
simplified  considerably  due  to  the  fact  that  =  0. 

Thus  Eq  (C-26)  has  the  following  terms  for  a  half-sin  load 
distribution: 


K .  .  =  & 
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-12 

6h 

6h 

4h2 

-6h 

2h 

-12 

-6h 

12 

-6h 

6h 

2h2 

-6h 

4h 

.1 


F. 

i 


PAE 

IT 


h 


< 


(CL  +  CO) 

+ 

frh 

(SL-SOP 

(CL+2C0) 

+ 

6E 

n 

(SL-SO) 

(CL  +  CO) 

- 
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+ 
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(SL-SO)^ 

> 


where 


CL  =  cos  —  ,  CO  =  cos  -  , 


E 


SL  =  sin  ~  ,  SO  =  sin 

E  E 


Sample  Problem  C-3.  A  uniform  clamped-free  rod  is 
subjected  to  a  half -sin  distribution  of  uniaxial  tension. 

A 

This  problem  is  discretized  into  two  elements  as  illustrated 
in  Fig  C-l.  The  reduced  global  equations  as  formed  by  the 
least  squares  technique  are 
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4h2  -6h  2h5 
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The  solution  to  these  equations  is 


•  117 


0 
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u2 
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U3 
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.2605 

U4 

f  AE  \ 

.3183 

U5 
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which  agrees  with  the  exact  solution. 


Appendix  D 

Derivation  of  Finite  Element  Equations 
for  the  Beam  in  Bending 


A  beam  in  bending  has  a  fourth-order  linear  differ¬ 
ential  equation  governing  its  transverse  displacement.  The 
equation  is 


AL  jI(x)E  aioixlL  P(x) 

dx2  dx2 


(D-l) 


If  w(x)  is  approximated  by  w(x)=Nj(x)W^  ;  where  j=l,2,...,M; 
the  residual  error  in  the  discretized  domain  will  be  of  the 
form 


a2  (  a -n. 

e  (x)  =  -  ]l(x)E  - w.  -  F(x)  #  0 

dx2  I  dx2  >  3 

i=l,2, . . . ,M  (D-2) 

The  method  of  weighted  residuals  is  also  used  for  this  prob¬ 
lem  to  derive  the  finite  element  equations.  The  integral 
equations  associated  with  each  of  the  following  finite  ele¬ 
ment  techniques  are  repeated  here  for  convenience. 
Collocation: 

e (x) 6 (x-x^)dx  =  0  (C-3) 
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Galerkin 


I 


e (x)N^ (x) dx 


0 


Least  Squares: 


I 


9e  (x) 
3w. 


1 


dx 


0 


(C-4) 


(C-5) 


Collocation 

Substituting  Eq  (D-2)  into  Eq  (C-3)  yields 


I 


dx4 


+  2 


dl  (x) 
dx 


d3N  . 
_ 1 

dx3 


+ 


dx2  dx2  J 


w 


j 


-  F  (x)  |  6(x-xi)dx  =  0 

j=l,2,...,M  (D— 3) 

The  transformation  from  a  global  coordinate  system  to  a 
local  coordinate  system  is  identical  to  that  for  the  axial 
rod  as  detailed  in  Appendix  C.  In  this  study  the  moment 
of  inertia  I(x)  was  assumed  to  be  either  constant  or  linear. 
That  is: 

I(x)  =  1(1- a  y*)  (0<a<1.0)  (D-4) 


Therefore 


and 


ft21  (-x)-  =  0  for  all  x 


dx 


Thus  the  integrated  Ie  associated  with  element  e  becomes 


'll1  - 1 


[  (e-1) h+x] /IE 
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0t  -rm 
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d3N . 
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i  w.-F(x)' 


=  0 


x=x. 


1*1/  2|  •  *  •  /in 
3“1 / 2 / • ■ • ,m 
e=*l,2 , . . .  ,E  (D-5) 


where  the  local  shape  functions  must  now  be  such  that 


d4N . 


dx" 


^  is  non-zero.  The  functions  chosen  are  the  septic 


shape  functions  (m=8)  for  a  C3  line  element.  (See  Appendix 
E.)  The  elemental  equation  for  element  e  can  be  written  in 
the  form 


I  =  K .  .w .  -  F.  =  0  (D-6) 

e  13  ]  1 

where  K.^  is  defined  for  section  as  ;  K^*4* 

is  the  part  of  the  stiffness  matrix  associated  with  the 
d“N. 


-i  terms;  and  K . . ^ 3  ^  is  the  part  associated  with  the 
dx*  13 

d3N. 

— — 1  terms. 
dx3 

Let  Ii  “  {l  “  Z  I  (e-Dh+^1  j-  and  x  =  ^  ,  then: 
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14 
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2axE 
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17 
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18 
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1-360  +  2340x  - 
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[  -4  +  3 Ox  - 
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?80x2  -  420x3 ] 

60x2  +  35x3 ] 


The  force  vector  F^  is  given  by  =  P  (constant) , 
i=l,2,...,8,  for  the  uniform  load  distribution,  and 
F^  =  P  sin  j  (e-l)h  +  x^|  for  the  half-sin  load  distribu¬ 
tion.  As  with  the  axial  rod,  the  collocation  points  x^ 
have  been  chosen  as  the  equally  spaced  points;  x^  *  -^y^-h, 

l*—  1 ,2,...  ,8. 


Galerkin 

For  the  Galerkin  method  Eq  (D-2)  is  substituted 
into  Eq  (C-4)  and  the  stiffness  term  is  integrated  by  parts 
twice.  The  result  is; 
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j-1,2, . . . ,M  (D-9) 


The  two  boundary  terms 


El  (x) 


d2w 

dx2 
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and 


EI(X)  +  E  smi.  sh> 

, 3  OX  j„2 


dx3 


dx 


L 

0 


both  vanish  for  the  boundary  conditions  studied.  The  fol¬ 
lowing  conditions  give  this  result.  Since  the  cubic  shape 

functions  satisfy  the  clamped  end  geometric  boundary 
clw 

conditions  (w  ~  ^  =  0) »  the  boundary  terms  vanish  for 
clamped  ends.  The  boundary  conditions  for  a  free  right 
end  are 
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d2  w 
dx2 
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I  (x)E 


d3  w 
dx3 


x=L 


0 


These  conditions  correspond  to  zero  moment  and  zero  shear 
force  respectively.  Therefore  if  the  right  end  is  free 
the  boundary  terms  again  vanish.  Equation  (D-9)  is  thus 
simplified;  and  if  the  same  coordinate  transformation  is 
performed  as  for  the  axial  rod,  the  result  for  element  e 
can  be  written  in  the  form  of  Eq  (D-6)  where 
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(D-10) 


and  for  the  uniform  load  distribution  F(x)  =  P  (constant) 
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(D-ll) 
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The  half-sin  load  distribution  F(x)  =  P  sin  (^)  gives  an 

Xj 


equivalent  force  vector  of 
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tt  (e-1) 


Least  Squares 

For  the  least  squares  method  Eq  (D-2)  is  substituted 
into  Eq  (C-5)  and  the  result  is: 
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d2I (x)  _ 

dx2 


Since  I(x)  *  1(1-  a |) 
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0 


If  Eq  (D-13)  is  written  in  local  coordinates  for  a  domain 
that  is  discretized  into  E  elements,  the  integral  for  ele¬ 
ment  e  is 


die4 


+  E 


+  E 


dl(ie) 

dx 

dl  (ie) 
dx 


F(x) 


EI(ie) 


0 


.  —  1 , 2  ,  • « • ,  m 
3-1,2, « •  • , m 
e=l ,2 (D-14) 


For  this  method  the  septic  shape  functions  are  chosen  (see 
Appendix  E) ,  and  let 


d  N.  x  x2  x3 

-g?  -  «i  +  tRi  +  !  Bi  +  T  Ci  *  V  Di’ 


Ai  +  Bix  +  C^x2  +  Dix3 


(D-15) 


where 


=  [0/  0/  0,  X/  0/  0^  0/  0] 

A.  =  —  [-840,  -480h,  -120h2 ,  -16h3,  840,  -360h,  60h2,  -4h3l 
1  h4 

B.  =—  [10080,  5400h,  1200h2 ,  120h3,  -10080,  4680h,  -840h2,  60h3] 
1  h4 


'H 


c.-—  [-25200,  -12960h,  -2700h2,  -240h3,  25200,  -12240h,  2340h2,  -180h3] 


1  h4 


D.  =—  [16800,  8400h,  1680h2,  140h3,  -16800,  8400h,  -1680h2,  140h3] 
1  h4 


Let  I(i)=  IL(1  -  |)  +  IR(f) 


Then  Eq  (D-14)  can  be  written  in  the  form  of  Eq  (D-6) 
where 


Kij  ■  ®>  {8<1L-2ILIR+IR,“it'; 


+  5 


+  ?  (ILZ-  6ILV5IR2|,C‘iCj+Ci“j> 


+  TS  'W1.  "Wi1 


+  (IL  -3ILrR+3IR  )AiAj 


+  *15  IL2"  5  ILIR  +  5  IR  )BiBj 
+  IL!"  2l  ILIR+6J  XR  ,CiCj 


+  (—  I  2-  —  I  I  +—  I  2)D.D. 
+  '56  XL  8  XLXR  +  4  XR  '  i  j 


+  IL2-  ILIR  +  |  IR2) 


+  (|  IL2-  |  ILIR+IR2)  (A^j+C^j) 
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+  'nV-l'A+IVVAV 


♦'iV-B  Vr-I  V»<BicfiBj' 


+  <ra>  V‘  M  Vr  + 1  V'  (BiDj+DiBj' 


+  *112  V-  168  It,IR+  16  IRI 1  (CiDj+DiCj*} 


3  T  2  ^2?  T  T  ,  5  T  2 


i—l r2r...,8 

j— 1 f 2» • « • |8 

e=l ,2, ...  ,E  (D-16 ) 


The  force  vector  for  the  uniform  load  case  F (x)  =  P  is  given 
by 

Ft  =■  Ph  f«R-ILl2ai  ♦  (3Ir-Il> 


B.  C.  D. 

+  (4Ir-Il)  -f  +  (5Id-It)  +  (6Id-It)  ^ 


R  L  12 


LR  L  20 


1—1 r  2  r •  •  •  r  8 

e=l,2 , . . . ,E  (D-17) 

The  uniform  beam  has  IL=IR=I  and  a=0;  therefore  Eq  (D-16) 
reduces ,  to 


Kij  -  EIShjAiBj+i  +i  (BiB.+AiCj+CiAJ) 

+  T  <BiCj4CiBj+AiDj+DiAj)  (CjC ^B^+D^ ) 


+  j  (C.D.+D.C.)  D,D. 
6  1313  713 


1—1  R  2  R  •  •  •  R  8 

3=1 r2r • • . r8 

e=l,2,...,E  (D-18) 
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and  Eq  (D-17)  reduces  to 


I- 


pi  -  PIh  {Ai  +  I  Bi  +  I  Ci  +  ?  Di}  (D-19> 

Substituting  the  appropriate  constants  in  Egs  (D-18)  and 
(D-19)  gives  the  result  for  Ig  shown  in  Eq  (D-20) .  The 
half-sin  load  distribution  results  in  an  equivalent  force 
vector  which  was  quite  unexpected.  The  force  vector  is 
given  by  the  following  equation 


F.  =  K. .E . 
i  ID  D 


i=l,2, ... ,8 

3=1 / 2  f • • . , 8 


where 


=  uniform  beam  stiffness  matrix  from  Eq  (D-20) 


and 


E  .  =  — 
D  * 


¥ 

sin 

ir(e-l) 

E 

¥ 

cos 

ir(e-l) 

E 

_  h 

~  ir 

sin 

ir(e-l) 

E 

cos 

ir  (e-1) 

E 

¥ 

sin 

it  e 

E 

¥ 

cos 

Tre 

E 

_  L 

*rr 

sin 

•ne 

II 

E 

cos 

tre 
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(D-20) 


Appendix  E 

Shape  Functions  and  Several  Derivatives 

The  shape  functions  used  for  a  particular  finite 
element  technique  determine  the  nature  of  the  approximate 
solution  (Ref  6:131) .  These  functions  as  used  in  this 
study  can  be  formed  from  the  familiar  Hermite  Polynomials 
(Ref  6:152-155).  Unfortunately,  this  method  does  not  con¬ 
vey  the  physical  nature  of  the  functions. 

A  linear  approximation  of  the  field  variable 

-  tal/ 

within  an  element  is  given  by  4>  (x)  =  [1  x]  <  >.  For 

a  line  element,  as  depicted  in  Fig  E-l,  the  nodal  parameters 
(j)^  and  <J>2  are  given  by 

i:;l  ■  c 


*2 

_ x 

2 

Fig  E-l.  C°  Line  Element 
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Inverting  this  relation  and  substituting  a^  and  aj  into 
the  approximation  equation  gives 


** 

1 

0 

*1  J 

4><x)  =  [1  x] 

1 

1 

•*2  j 

“  h 

h 

The  product  of  the  two  matrices  is  defined  as  the  vector  of 
shape  functions  N ^ ,  where  j  denotes  the  shape  function 
associated  with  nodal  parameter  The  approximation  equa¬ 

tion  is  now  of  the  form: 


~  m 

<Mx)  = 

j=l 


j — 1  /  2  f  •  •  .  ,m 


where  m  =  number  of  nodal  parameters  associated  with 
the  element. 

This  type  of  element  is  said  to  have  C°  continuity  because 
only  the  function  is  continuous  between  two  connected  ele¬ 
ments.  The  linear  shape  functions  are  given  in  Fig  E-2. 
These  shape  functions  exhibit  properties  of  the  "kronecker" 
delta  function.  That  is: 


N 


j 


i=j 


A  line  element  possessing  continuity  will  have 
continuous  first  derivatives  of  the  field  variable  between 
two  connected  elements.  This  requires  a  cubic  approxima¬ 
tion  of  the  field  variable.  By  duplicating  the  process 
followed  for  the  linear  shape  functions  one  can  transform 
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Fig.  E-2 .  Linear  Shape  Functions  (Ref  5:154) 
the  equation  <Mx)  =  +  a2*  +  a3**  +  a4^3  into 


4 

<P  (x)  =  ^Nj^j 

j=l 


Fig.  E-3.  Line  Element 


The  cubic  shape  functions  associated  with  the  line  ele¬ 
ment  are  given  in  Fig  E-4.  The  linear  and  cubic  shape  func 
tions  are  summarized  in  Table  E-I. 


Fig.  E-4.  Cubic  Shape  Functions  (Ref  6:154) 

The  highest  order  element  used  in  this  thesis  is  a 
C3  line  element.  It  allows  continuity  of  the  third  deriva¬ 
tive  of  the  field  variable  between  elements.  The  approxi¬ 
mation  equation  is  $  (x)  =  +  a2*  +  a3^2  +  a4^3  +  a5*1’ 

+  agX5  +  a^x6  +  a8x7.  Since  there  are  eight  undetermined 
parameters  the  C3  line  element  must  be  defined  with  four 
nodal  parameters  at  each  of  its  two  nodes.  Thus  the 


vector  is 


in 


Table  E-l 

Linear  and  Cubic  Shape  Functions  with 
Derivatives  of  Specific  Interest 
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*.--W 


The  third  and  seventh  parameters  are  proportional  to  the 
bending  moment  in  a  beam  at  x=0  and  x=h.  .  The  fourth  and 
eighth  parameters  are  proportional  to  the  transverse  shear 
force.  The  septic  shape  functions  are  given  in  Table  E-II, 
along  with  their  fourth  derivatives.  Shape  functions 
through  N4  are  shown  in  Fig  E-5;  the  other  four  shape  func¬ 
tions  are  simply  a  reflection  about  the  point  x=h/2. 


Fig  E-5.  Septic  Shape  Functions  (Ref  5:155) 
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Table  E-II 


Septic  Shape  Functions  with  Derivatives 
of  Specific  Interest 


l-35x4 

+ 

85xS 

- 

7  Ox6 

+ 

20x7] 

h[x-20x4 

+ 

45xs 

- 

36x6 

+ 

10x7  ] 

h2  [  |x2-  5x4 

+ 

10x5 

- 

i^x6 
2  X 

+ 

2x7] 

h3  l|x3-|x4 

+ 

x5 

- 

2  6 
3  X 

+ 

I1*7 ) 

35x4 

- 

84x5 

+ 

70x6 

- 

20x7  ] 

h[-15x4 

+ 

39x5 

- 

34x6 

+ 

10x7] 

h=[|x* 

- 

7x5 

+ 

13xs 
2  X 

- 

2x7  ] 

+ 

lx5 
2  x 

- 

Ixs 
2  X 

+ 

?*7i 

-[-840  +  10080x  -  25200x2  +  16800x3] 
h4 

-[-480  +  5400x  -  12960X2  +  8400x3] 

h3 

-[-120  +  1200x  -  2700x2  +  1680x3] 

h2 

-jj-l  -16  +  120x  -  240x2  +  140x3] 

— [  840  -  10080x  +  25200X2  -  16800x3] 

h4 

-[-360  +  4680x  -  12240X2  +  8400x3] 

h3 

— [  60  -  840x+  2340x2  -  1680x3] 

h2 

-4  +  6 Ox  -  180x 2  +  140x  *] 


Appendix  F 

Derivation  of  Finite  Element  Equations  for  the 

'  t 

h  i  Nonlinear  Problem  with  Sample  Calculations 

The  differential  equation  given  for  the  nonlinear 
problem  studied  is 

U-^xx  =  0  0<x<l  (F-l) 

with  the  boundary  conditions 

<M0)  =  0 
*x(0)  =  1 

4>(1)  =  (constant)  (F-2) 

where  the  subscript  x  implies  differention  with  respect  to 
the  independent  variable  x.  Equation  (F-l)  can  be  stated 
in  two  different  but  equivalent  forms.  These  are 

A»  *x  -  1  '*x2lx  -  0  (r-3> 

B)  [(l-*x>2]x  =  0  (F-4) 

The  domain  was  discretized  into  two  elements  of 
length  h^  and  ^respectively,  with  the  global  nodal  param¬ 
eters  given  in  Eq  (F-5) .  The  two  elements  are  joined  at  a 
double  node  as  shown  in  Fig.  F-l. 
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J 


=  <M0) 


*  =  d<(>  (0) 

v2  d  x 


<fi  3  =  *0^  ) 


♦4  - 


d4>(h1  ) 


(h1+) 
d<J>  (h1+) 


=  <f>  (h1+h2) 


<3-<f)  (1) 
dx 


(F-5) 


where  $3  and  are  the  approximate  solutions  to  the  dis¬ 
placement  (potential  function)  and  the  slope  (velocity) 
respectively  at  the  right  end  of  the  first  element;  and 
<j>5  and  <J>g  are  the  approximate  solutions  at  the  left  end  of 
the  second  element.  The  discretized  domain  is  shown  in 


Fig.  F-l. 


Fig  F-l.  Discretized  Domain  for  Nonlinear  Problem 
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It  is  necessary  to  have  at  least  one  double  node  because 
the  differential  equation  with  the  boundary  conditions  given 
by  Eq  (F-2)  has  exactly  one  discontinuity.  The  double  node 
allows  the  specification  of  the  discontinuity,  or  jump,  in 
the  approximate  solution.  The  jump  is  modeled  by  two  equa¬ 
tions  .  They  are : 


*3  =  *5 

♦4  =  +  A  (F-6) 

where  A  is  the  change  in  slope  between  elements.  The  first 
equation  insures  continuity  of  the  function  across  the  inter¬ 
element  boundary;  and  the  second  guarantees  the  required 
change  in  slope. 

With  this  discretization,  the  finite  element  equa¬ 
tions  can  be  derived  by  the  same  process  as  the  linear  prob¬ 
lems.  (See  Appendices  C  and  D.) 

Equation  (F-6)  is  solved  hy  an  iterative  process. 

With  the  substitution  <J>n(x)  =  Nj(x)<|>jn  ,  the  iterative 
equation  for  element  e  becomes 


=  0 


0<x<h  i=l,2,3,4 

e=l,2  j=l,2,3,4 

k=l ,2, 3,4  (F-7) 


« 

i 


\ 
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where  <j>jn  are  the  nodal  solutions  obtained  from  the  n — 
iteration. 

By  the  exact  same  process  Eq  (P-4)  gives 


/•h 

N.  N  .  dx<t>  .n+1  =  / 

Xx  ^  3  J0 


(l-Nk_4>kn)Ni_d^ 


0<x<h  i=l,2,3,4 

e=l/2  j=l,2,3,4 

k=l ,2 ,3,4  (F-8) 


The  shape  functions  chosen  are  the  cubic  functions 
given  in  Appendix  E.  The  local  shape  function  derivatives 
can  be  written  in  the  form 

0<x<h  i=l ,2,3,4  (F-9) 


where 


Ai  “  >-E’  *4'  5'  -21 


Bi  "  [  h'  3'  '  h'  31 ‘ 


With  N._  in  this  form  Eqs  (F-7)  and  (F-8)  can  be 
1x 

written  in  the  following  form  for  h^  =  hj  =  h 


Kij^jn+1  =  Fi 


isl $  2  f  •  • .  ,4 
j=l f 2| ■ • . • 4 
e=l  ,2 


where 


(F-9) 


Kij  -  410  [AiAjr1+(AiBj+BJAj)r2  +  B.B^J 

rl  “  84&1  +  105e2  +  140 
T 2  ~  7081  +  84b2  +.105 


r3  =  eOBj^  +  7062  + 


84 


and  for  Eq  (F-7)  -  formulation  A 

Si  "'A  I2l*1%3n)  +»  42%4n>) 


S2  -  -jj-  C3 +h  (2+2%4n)] 


FA  =  0  i-1,2,. 

and  for  Eq  (F-8)  -  formulation  B 

@1  =-h  t2(4>!%3n)  +h{^2n+(j,4n)] 

82  -  \  [3  (<J>1rl-<t)3n)  +h (2$2n+4>4n)  ] 

Fi  (18^i  +  3082  +  60) 

F2  «  (2481  +  3562  +  60) 

'3  *  To  (1861  +  3062  +  60) 

P4  '  S5  ,681  +  5B2» 

The  products  of  A^  and  in  matrix  form  are  as 
h  -  1/2. 


(F-10) 


,4  (F-ll) 


(F-12) 


follows  for 
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V 

144  48 

-144 

24 

[AAl  =  [A^J  = 

16 

-48 

8 

144 

-24 

(symmetric) 

4 

-288  -84 

288 

-60 

-24 

84 

-18 

[AB]  -  [A.B.+B.A.]  = 

l  j  l  j 

-288 

60 

(symmetric) 

12 

144  36 

-144 

36 

[BB]  =  IB±B^]  = 

9 

-36 

9 

144 

-36 

(symmetric) 

9 

Sample  Problem  F-l.  The  nonlinear  problem  for 
A=2  is  solved  by  formulation  A.  The  initial  guess  is 
♦?  =  [0,  1,  1,  1 ,  1,  -1,  0,  -11 

For  element  1 : 


K. 


*i  - 

3, $2  =  -3 

r  i  - 

77, r2  =  63, 

r3  =  54 

thus, 

720  348 

-720 

12 

1 

206 

-348 

-32 

840 

720 

-12 

(symmetric) 

38 

143 


For  element  2: 

($2  ~  ~ 3 ,  B 2  ~  3 

rx  =  203,  r2  =  i 


=  3 

=  147 

'  r3  * 

114 

1500 

-3312 

156 

746 

-1500 

4 

3312 

-156 

.ric) 

74 

Assembling  the  elemental  matrices  in  the  global  system  equa 
tions  gives 


748 

-720 

12 

206 

-348 

-32 

720 

-12 

38 

1500 

-3312 

156 

746 

-1500 

4 

3312 

-156 

74 

4} 


[symmetric) 


Substituting  the  jump  conditions  and  the  boundary  condi¬ 
tions  gives  the  reduced  system  equations.  They  are 


720  -12  0  0  0 

-12  38  0  0  0 

1  0-100 

0  10-10 

0  0  156  4  74 


r*n 

348 

*4 

32 

*5 

►  =  < 

0 

*6 

2 

1*8  j 

0 

Solving  this  equation  gives  the  following  which  matches  the 
exact  solution. 

♦£  -  [0,  1,  0.5/  1,  0.5,  -1,  0,  -1] 

The  next  iteration  is  performed  with  <)Kn  =  <J>?  which  gives 
.  This  second  iteration  gives  the  same  numerical 
results  as  the  first,  implying  convergence  of  the  algorithm. 


Collocation 

The  domain  is  discretized  in  the  same  manner  as 
for  the  Galerkin's  method.  The  derivation  of  the  finite 
element  equations  is  by  the  same  procedure  as  followed  in 
Appendices  C  and  D  for  this  method.  The  two  formulations 
give  the  following  global  equations  for  element  e. 

Formulation  A: 


_ 

x=x. 


n+1 


0 


0£x<h  i=l,2,...,m 

e=l,2  3&lf2,...,m 

k=l,2, . . . ,m  (F— 13) 
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Formulation  B 


u-V-V>Hi,,+VV~An 


'xx  J  x  xx 


4,  n+1 

_  _  j 

x=x^ 


x=x. 


0£x£h  i=l,2,...,m 

e=l, 2  j=l  ,2 , . . . ,m 

k=l,2 , . . .  ,m  (F-14) 


The  cubic  shape  functions  given  in  Eq  (F-9)  are  used  for 
this  method  also. 
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